#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
#include <omp.h>

int criterio_linha(double **,int);

int main(int argc,char *argv[]){

	int J_ORDER, J_ITE_MAX, J_ROW_TEST;  /**Variaveis do programa**/
	double J_ERROR;
		
	int sizeBlock, finalize = 0, numRank, numProc, nameLen;
	int count = 0,i,j,aux = 0,stop = 0;
	double div,dif,result;
	double numerator = 0.0, denominator = 0.0;

    double **matrixA;
	double *VB;
	double *X;
	double *X1;
	double *MA;
	double *buffer;	

	char nameProc[MPI_MAX_PROCESSOR_NAME];

	//MPI_Status status;
    MPI_Init(&argc, &argv);
	/**Retorna o rank do processo**/
    MPI_Comm_rank(MPI_COMM_WORLD, &numRank);
	/**Retorna a quantidade de processos**/ 
    MPI_Comm_size(MPI_COMM_WORLD, &numProc);
	/**Retorna o nome do processador**/
	MPI_Get_processor_name(nameProc, &nameLen);

	printf("Processor's Name: %s | Rank: %d\n", nameProc, numRank);

	if(numRank == 0){/**O processador de rank 0 ira parsear o arquivo de entrada, pegando os dados coletados
						e distribuindo nas matrizes e variÃ¡veis correspondentes **/

		scanf("%d",&J_ORDER);    /**Ordem da matriz**/
		
		if (J_ORDER % numProc != 0){/**Somente aloca os processos se conseguir dividir as matrizes entre eles**/
       		printf("Nao foi possivel dividir a matriz entre os processos alocados.\n");
       		finalize = 1;/**Varivel que indica se os processos tem que ser finalizados**/
    	}		

		scanf("%d",&J_ROW_TEST); /**Local a ser realizado o teste da linha**/
		scanf("%lf",&J_ERROR);   /**Erro permitido**/
		scanf("%d",&J_ITE_MAX);  /**Numero maximo de iteracoes**/			
	
		MA = (double *)malloc(sizeof(double)*J_ORDER*J_ORDER);/**Aloca o array MA[]**/		
		
		matrixA = (double **)malloc(sizeof(double *)*J_ORDER); /**Aloca a matriz matrixA[][], usada no teste de linha**/
		for(i = 0;i < J_ORDER;i++)
			matrixA[i] = (double *)malloc(sizeof(double)*J_ORDER);

		VB = (double *)malloc(sizeof(double)*J_ORDER);    /**Aloca o vetor VB[]**/

		for(i = 0;i < J_ORDER;i++){  /**Coloca os dados do arquivo na matriz matrixA e no array MA**/
			for(j = 0;j < J_ORDER;j++){
				scanf("%lf",&matrixA[i][j]);
				MA[aux] = matrixA[i][j];
				aux++; 
			}
		}			

		for(i = 0;i < J_ORDER;i++)/**Coloca os dados do arquivo no vetor VB**/
			scanf("%lf",&VB[i]);
		
		X = (double *)malloc(sizeof(double)*J_ORDER);/**Alocando os vetores soluÃ§Ãµes X[], X1[] e o vetor auxiliar buffer**/
		X1 = (double *)malloc(sizeof(double)*J_ORDER);
		buffer = (double *)malloc(sizeof(double)*J_ORDER);

		if(!criterio_linha(matrixA,J_ORDER)){/**Se a matriz nao passar pelo criterio da linha, finaliza os processos**/
			printf("NÃ£o obedece ao criterio da linha.\n");
			finalize = 1;
		}
		
		sizeBlock = J_ORDER / numProc;/**Calcula o tamanho do bloco da matriz que cada arquivo vai computar**/
	
	}

	MPI_Bcast(&sizeBlock, 1, MPI_INT, 0, MPI_COMM_WORLD);/**Os dados que serÃ£o usados pelos nÃ³s sÃ£o compartilhados aqui**/
	MPI_Bcast(&finalize,1,MPI_INT,0,MPI_COMM_WORLD);
	MPI_Bcast(&J_ORDER, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&J_ERROR, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
	MPI_Bcast(&J_ITE_MAX, 1, MPI_INT, 0, MPI_COMM_WORLD);

	if(finalize == 1){/**Finaliza os processos**/
		MPI_Finalize();
		exit(-1);
	}

	if(numRank != 0){

		MA = (double *)malloc(sizeof(double)*sizeBlock*J_ORDER); /**Aloca a matriz MA[][], o vetor VB[], e os vetores solucoes X[] e X1[]**/
		VB = (double *)malloc(sizeof(double)*sizeBlock);    
		X = (double *)malloc(sizeof(double)*J_ORDER);
		X1 = (double *)malloc(sizeof(double)*sizeBlock);
	
	}
 	
	for(i = 0;i < J_ORDER;i++)/**Inicializa o vetor X[]**/
			X[i] = 0.0;

	MPI_Scatter(MA, sizeBlock * J_ORDER, MPI_DOUBLE, MA, sizeBlock * J_ORDER, MPI_DOUBLE, 0, MPI_COMM_WORLD);/**Divide o vetor MA[] e o vetor VB[] para cada no**/
    MPI_Scatter(VB, sizeBlock, MPI_DOUBLE, VB, sizeBlock, MPI_DOUBLE, 0, MPI_COMM_WORLD);/**jÃ¡ repartindo os dados**/

	for(count = 0; stop == 0 && count < J_ITE_MAX;count++){/**FunÃ§Ã£o de Gauss-Jacobi**/
		denominator = 0.0;
		numerator = 0.0;

		omp_set_num_threads(4);
		#pragma omp parallel for default(none) shared(X1, MA, J_ORDER, X, sizeBlock, numRank, VB) private(i, j)
		for(i = 0;i < sizeBlock;i++){
			X1[i] = 0.0;
			for(j = 0;j < J_ORDER;j++){
				if(i*J_ORDER + j != (J_ORDER + 1) * i + sizeBlock * numRank)
					X1[i] = X1[i] + MA[(i * J_ORDER)+j] * X[j];
			}
			X1[i] = 1 / MA[(J_ORDER + 1) * i + sizeBlock * numRank] * (VB[i] - X1[i]);
		}		

		if(numRank == 0){/**Copia o vetor X[] em buffer[] para o calculo do erro mais a frente**/
			for(i = 0;i < J_ORDER;i++)
				buffer[i] = X[i];
		}

		MPI_Allgather(X1, sizeBlock, MPI_DOUBLE, X,sizeBlock,MPI_DOUBLE, MPI_COMM_WORLD);/**Reune as solucoes encontradas por cada no no vetor X[]**/

		if(numRank == 0){/**Calculo do erro**/
			#pragma omp parallel for default(none) shared(J_ORDER, X, X1, buffer, numerator, denominator) private(i, dif)
			for(i = 0;i < J_ORDER;i++){
				dif = X[i] - buffer[i];

				if(dif < 0)
					dif *= (-1);
		
				#pragma omp critical
				if(numerator < dif)
					numerator = dif;

				dif = X1[i];

				if(dif < 0)
					dif *= (-1);

				#pragma omp critical
				if(denominator < dif)
					denominator = dif;		
			} 
			div = numerator/denominator;

			if(div < J_ERROR)
				stop = 1;
		}		

		MPI_Bcast(&stop,1,MPI_INT,0,MPI_COMM_WORLD);/**Reparte a informaÃ§Ã£o se o o nÃ³ deve parar o processamento ou nÃ£o**/
	}

	if(numRank == 0){/**Calculo e apresentacao dos resultados finais**/
		if (count < J_ITE_MAX){                      
        	for (i = 0; i < J_ORDER; i++)
            	result += MA[(J_ROW_TEST * J_ORDER) + i] * X[i];
            printf("Iterations: %d\n", count);
            printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, result, VB[J_ROW_TEST]);
        }
		else
			printf("Iteration overflow!\n");

		free(matrixA);/**Libera os vetores usados**/
		free(X);
		free(X1);
		free(MA);
		free(VB);
		free(buffer);
	}

	MPI_Finalize();

	return 0;
}

int criterio_linha(double **M,int J_ORDER){/**Calculo do critÃ©rio da linha, para ver se a matriz ira convergir**/
	int i,j;
	double sum = 0.0,max = 0.0;

	for(i = 0;i < J_ORDER;i++){
		sum = 0.0;
		for(j = 0;j < J_ORDER;j++){
			if(i != j)
				sum += M[i][j];
		}
		sum = sum/M[i][i];
		if(max < sum)
			max = sum;
	}

	if(max < 1)
		return 1;
	else
		return 0;
}
